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Abstract 


The decay of isotropic turbulence in a compressible flow is examined by direct nu- 
merical simulation (DNS). A priori analyses of the DNS data are then performed to 
evaluate three subgrid-scale (SGS) models for large-eddy simulation (LES): a gener- 
alized Smagorinsky model (Ml) [J. Fluid Mech., 238, 1 (1992)], a stress-similarity 
model (M2) [«/. Fluid Mech., 275, 83 (1994)], and a gradient model (M3) [Theoret. 
Comput. Fluid Dynamics, 8, 309 (1996)]. The models exploit one-parameter second- 
or fourth-order filters of Pade type, which permit the cutoff wavenumber k c to be 
tuned independently of the grid increment Ax. The modeled (M) and exact (E) SGS- 
stresses are compared component-wise by correlation coefficients of the form C(E, M) 
computed over the entire three-dimensional fields. In general, Ml correlates poorly 
against exact stresses (C < 0.2), M3 correlates moderately well (C ~ 0.6), and M2 
correlates remarkably well (0.8 < C < 1.0). Specifically, correlations C(E, M2) are 
high provided the grid and test filters are of the same order. Moreover, the highest 
correlations (C « 1.0) result whenever the grid and test filters are identical (in both 
order and cutoff). Finally, present results reveal the exact SGS stresses obtained by 
grid filters of differing orders to be only moderately well correlated. Thus, in LES the 
model should not be specified independently of the filter. 

Keywords: large-eddy simulation, subgrid-scale modeling, isotropic turbulence, compress- 
ible flow, digital filtering 
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1 Introduction 

By definition, direct numerical simulation (DNS) is the numerical solution of the Navier- 
Stokes equations without recourse to empirical models. In concept, the fluid motions are 
resolved down to the Kolmogorov wavenumber k v , at which eddies succumb to viscous dissi- 
pation. In general, computational workload for fully-res^^i bNS scales as Me 3 , where 
Re is the Reynolds number. For the complex, high-Reynolds-number flows of engineering in- 
terest, the computational requirements of DNS are staggering and prohibitive. Consequently, 
large-eddy simulation (LES) remains an attractive (and necessary) alternative. 


In contrast to DNS, LES is based on a decomposition of the flow-field variables into 
resolved and unresolved (subgrid) scales. The decomposition is effected by a spatial grid 
filte ring operation with an associated cutoff wavenumber k e . The grid-filtered Navier-Stokes 
equations incorporate the subgrid-scale (SOS) stress tensor t#, which quantifies the interac- 
tions between the smaller (k > k c ) and larger (k < k c ) scales of motion. The larger scales 
are resolved in space and time on a suitable computational grid; however, the SGS stresses 
are modeled or otherwise approximated. 

The relationship between DNS and (ideaiized) r LES is illustrated in Fig. I, which presents 
the Kolmogorov energy spectrum. The inertial range of the spectrum is characterized by 
energy decay at the rate E(k) ~ k -5 A. By implication, r w varies with the cutoff k c \ that is, 
T k i = r k i(k c ). For DNS, k c = k v and, for all practical purposes, r w (k„) = 0. 


A general criticism of LES [1} is that the truncation errors of the numerical method often 
contaminate the SGS dissipation. There is growing realization (e.g., Vasilyev et al. [2]) that 
this situation can be avoided by independent specification of the grid resolution A®* and the 
wavenumber cutoff k*. (Throughout this work, asterisks denote dimensional quantities;) In 
our view, k* should be specified on the basis of physical considerations; that is, k* c should lie 
in the inertial range of the Kolmogorov spectrum (Fig. 1). On the other hand, Ax* should be 
determined by numerical considerations; that is, by the grid resolution necessary to resolve 
the smallest eddies (those of wavelength A* = 2 tt / k*) for the numerical scheme of choice. 
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For example, spectral, sixth-order, fourth-order, and second-order difference approximations 
typically require 2, 10, 16, and 32 gridpoints, respectively, per wavelength of the smallest 
resolved scale. If we define the dimensionless wavenumber a = k* Ax*, then a c = Ax* 
defines the dimensionless filter cutoff. It is highly desirable, therefore, to exploit continuously 
tunable filters, whose cutoffs a c are not “hard-wired” to the grid resolution. For this purpose, 
we adapt one-parameter families of filters of Pade type, as described in Lele [3]. 

Our present work parallels the experimental investigations of Liu et al. [4] and the 
computational work of Vreman et al. [5]. Each paper evaluates the potential of candidate 
SGS models by means of a priori analyses. (The latter work also performs a posteriori 
analyses.) In their ground-breaking experiment, Liu et al. [4] investigated the structure of 
turbulence in the far field of an axisymmetric water jet by means of particle-displacement 
velocimetry. By planar filtering of the velocity data, they reconstructed (2D slices) of three 
components of the exact SGS-stress tensor and correlated various modeled stresses against 
their exact counterparts. To our knowledge, their experiment represents the first a priori 
analysis conducted with experimental rather than computational data. Two of the models 
considered by Liu et al. [4] were the standard (incompressible) Smagorinsky model and 
a stress-similarity model. In Vreman et al. [5], three candidate SGS models (and their 
dynamic c min t.ftrpa.rts) were evaluated by both a priori and a posteriori analyses, using 
data obtained by DNS of a weakly compressible turbulent mixing layer. Specifically, they 
examined Smagorinsky, similarity, and gradient [6, 7] models. 

In the present paper, we evaluate, for the case of decaying isotropic turbulence, the three 
models considered by Vreman et al. [5]. Our evaluation is limited to a priori analysis only 
and has a different focus than the work of Vreman et al. [5]. Whereas, they exploit a top-hat 
filter with a fixed filter width (twice the grid increment) and focus on the SGS models, we 
exploit Pade- type filters and focus specifically on the effects of the filter. In particular, we 
examine the effects of the cutoff and the order property of the filter on the fidelity of the 
models. 

The next section discusses our two candidate filters: second- and fourth-order filters of 
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Pade type with continuous variation of the cutoff a c . In Section 3, we present the unfiltered 
and filtered compressible Navier-Stokes equations, which serve as the governing systems for 
DNS and LES, respectively. Section 4 describes the three SGS models to be considered: 
the generalized Smagorinsky model of Erlebacher et al. [8], the stress-similarity model of 
Liu et al. [4], and the generalized Clark model of Vreman et al. [7]. Section 5 presents 
parameters of the test flow: decaying isotropic turbulence. The a priori analysis is presented 
in Section 6. In particular, we present correlations between the exact SGS-stress components 
and their modeled counterparts. A brief discussion of the results is presented in Section 7, 
and conclusions are offered in Section 8. 

2 One-Parameter, Low-Pass Filters 

Many types of spatial, linear, low-pass filters can be used for LES, the most popular being 
top-hat, Gaussian, and spectral filters. The properties of admissible filters can be found in 
many sources, including Erlebacher et al. [8]. 

Following Lele [3], we consider one-parameter families of filters of Pade type that are both 
symmetric and positive semidefinite. The positivity of these filters ensures that the SGS- 
stress tensor is realizable in the sense of Vreman et al. [9|. A filter is best understood in terms 
of its transfer function, which depicts its effect on modes of wavenumber a. The transfer 
functions of a family of second-order Pade filters are shown in Fig. 2 for selected values 
of the cutoff at c . The dissipation imposed by the filter varies continuously over the range 
0 < ot c < 7T, decreasing with increasing ot c . Specifically ot c = i r turns off the filter, a c = nj2 
results in the discrete top-hat filter (with weights at adjacent nodes of [1/4, 1/2, 1/4]), and 
at e « 0 results in extreme dissipation. In general, one-parameter filters are advantageous 
relative to fixed-width filters in that a c can bo adjusted so as to specify k c independently of 
Ax. For some purposes, we will also exploit a one-parameter, fourth-order Pade filter with 
the same parameter range. Details of the both filters are relegated to the Appendix. 

For our present purposes, multidimensional filtering is accomplished by filtering succes- 
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sively in each coordinate dimension. The same cutoff is exploited in each dimension. 


3 Governing Equations 


For DNS, the Navier-Stokes equations are solved numerically on a grid sufficiently fine to 
resolve all energetic scales of motion. 

_ . _ _ i / 

3.1 Compressible Navier-Stokes Equations 


In tensor notation, the dimensionless equations governing the evolution of an ideal gas in 
time t and Cartesian space x = [x\, 12,13] are 

dp djpUk) 
dt dxk 
djpUk) djpUkUi ) 
dt dxi 

d{Er) d\{Er + pW] 

dt dxk 

where 


dp t 1 do k i 
dxk Re di( 
dhk [ 1 d{ui<r kl ) 
dxk Re dxk 


Et — p{e + 


2 } 


e 


hk 


T 

7(7 — 1 )M 2 

, W dT 
(7 - 1 )M 2 PrRedxk 


( 2 ) 


and repeated indices imply summation. Here, p, p, T, e, Er, p, and k are density, ther- 
modjmamic pressure, absolute temperature, specific internal energy, total energy, dynamic 
viscosity, and thermal conductivity, respectively, and u = [111,112,1x3] is the velocity vec- 
tor. The quantities above have been rendered dimensionless by the choice of a refer- 
ence length scale L* and suitable reference values (denoted by subscript “r”) as follows: 
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(3) 


t = T’/r; 

« 

V 

II 
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From the nondimensionalizations above, four dimensionless parameters arise: the Reynolds 
number, the Mach number, the Prandtl number, and the ratio of specific heats. These are 
defined, respectively, as follows: ! 


Re •=* 
M = 
Pr = 
7 = 


Px 

«r_ 


vW 





( 4 ) 


Here, C* and C* are the specific heats at constant pressure and constant volume, respectively, 
and R* is the ideal gas constant. Throughout this work, we presume the constant values 
7 = 1.4 and Pr = 0.7. The viscous stress tensor &ki is defined 


<?ki = 2 fiSki - ~v.D8u 


( 5 ) 


where £*< is the Kronecker delta, D — ^ is the compressible dilatation, and Ski is the 
symmetric strain-rate tensor defined as 


Ski 


1 fduk . du t \ 

2 y dxi dx k ) 


( 6 ) 


Equations 1 describe the conservation of mass, momentum, and energy, respectively. The 
governing system is closed by imposing an equation of state for an ideal gas, namely 

pT 


P = ^ = <7-lVe 


(7) 


The variation of p with respect to T is modeled by Sutherland’s law, and k = p/Pr. 


3.2 Favre-Filtered Navier-Stokes Equations 


For LES of compressible flow, density-weighted (Favre) filtering arises naturally from the 
use of conservative variables. For example, the Favre-filtered velocity components, denoted 
by tildes, are given by 


u k = 


Wk 

p 

(pu k = piLk) 


( 8 ) 


The Favre-filtered continuity equation is form invariant; that is, 

dp d(pUk) 
dt dx k 


0 ) 


as is the filtered equation of state p = pf/^M 2 ). However, the momentum equation acquires 
the following filtered form: 


d(pUk) d(pUkUi) _ d p 1 ddki dr ki 

dt dxi dxk -Re dxi dx / 


( 10 ) 


where r k i is the SGS-stress (or the residual-stress) tensor. For compressible flow, the residual- 
stress tensor assumes the form 

T k l = P{u k Ul - UkUl) (11) 

As written, Eq. 10 is exact; however, in practice it becomes inexact because &ki, being 
unavailable, is usually approximated from Ski and ft, with residual stresses due to their 
nonlinear interactions being neglected. Following Erlebacher et al. [ 8 ], an alternate form of 
the dimensionless grid-filtered energy equation is 


g + (7-1*0 + 


dt ' dxk 
where the dissipation function is 


1 .J d 
M 2 PrRe dx k 




+ 2Z±$ + 

Re 7 M 2 dx k 


( 12 ) 


* = 2liSu ^ ‘ I" 02 


(13) 


and the residual heat-flux vector is 


Qk = p{u k T — UkT) 


(14) 
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Equation 12 is inexact because it has been derived by neglecting certain pressure-velocity 
terms, and because, in practice, the term is typically approximated from the resolved 
scales. 

4 Subgrid-Scale Stress (SGS) Models 

We consider three residual-stress models, denoted Ml, M2, and M3, as described, respec- 
tively, below. 

4.1 Generalized Smagorinsky Model (Ml) 

A generalization of the Smagorinsky model to compressibleflow is presented in Erlebacher 
et al. [8], which we adapt as follows. Let MRu denote the modeled residual-stress tensor 
such that 

MRki — (15) 

where §ij is the resolved-scale strain-rate tensor, and f*r is the eddy viscosity defined as 

W = 2C*A s 0|S| ; |S| = s/SuSu (16) 

Here, Cr is a model (instant, and A is a characteristic length scale related to the grid spacing. 
Typically A = \/Az 2 + Ay 2 + Az 2 . Similarly, the modeled residual heat-flux vector MQk 
is given by 

(17) 

where C q is a model constant related to Cr through the turbulent Prandtl number Prr- 

4.2 Scale-Similarity Model (M2) 

On the basis of their experimental observations, Liu et al. [4] propose (for incompressible 
flow) the simple stress-similarity model 

Tfei « ClCm (18) 
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where c L is a constant, and Cm is the resolved-turbulent-stress tensor. A computable quan- 
tity, Cm is extracted from the resolved fields by applying a second filter, termed the test filter 
and denoted by a hat. Specifically, 

C M = t k %,i - (19) 

In general, the cutoffs associated with the grid and test filters may be different. Bardina’s 
[10] model is the special case of Eq. 18 for which ci — 1 and r = 1, where r = A/A, and A 
and A are the filter widths associated with the test and grid filters, respectively. 

For compressible flows, attempts to distinguish Favre grid- and test-filtered quantities 
result in cumbersome notation. Therefore, we adopt the simplifying convention that, for any 
twice-filtered quantity, the first and second filtering operations are presumed to correspond 
to the test and grid filters, respectively, which may have different cutoffs. Following this 
convention, the resolved-turbulent-stress tensor is 

Cm = - u k ui) ( 20 ) 

and its thermal counterpart is 

q k = p(Z k f - fij) (21) 

4.3 Generalized Clark Model (M3) 

Using Taylor-series expansions of the filtered velocity fields in terms of the unfiltered veloc- 
ities, Vreman et al. [7] derive the gradient model 

'‘■ = -i^ljf| +0(A4) (22) 

(The negative sign is included for consistency with our definition of the residual stress 
(Eq. 11)). Their model generalizes that of Clark et al. [6] to compressible flows. Crit- 
icism of the use of Taylor expansions by Love [11] has been addressed in Vreman et al. 
[7]. Here, our purpose is not to defend or refute the Taylor-series approach, but simply to 
evaluate the proposed model by a priori analysis. 
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5 Test Case: Decaying Isotropic Turbulence 

We consider the decay of an isotropic turbulent compressible flow. Specifically, the test case 
corresponds approximately to case ife96 of Blaisdell et al. [12], whose parameter values are 
restated below for completeness. 

Re 
M 

T* 

-*00 

where 

Pv — Po 

= C$ = y/rirrs (24) 

r r * = r 0 * 

and the subscript “0” refers to spatially averaged values at the initial time. For reasons of 
computational efficiency, the computational domain is a cube of edge length 2#, with periodic 
boundary conditions in each coordinate dimension. The flow is initialized with the top-hat 
energy spectrum shown by solid line in Fig. 3, for which the phases are prescribed randomly. 
Whereas Blaisdell et al. [12] exploits a power law for the viscosity model, as mentioned 
previously, we employ Sutherland’s law. As a consequence, some of our initial parameters 
differ marginally from those published in Blaisdell et al. [12]. The initial turbulent Mach 
number M t = \j< >/\j'y < p > / < p > = 0.7. Here, angle brackets denote volume 

averages over the entire domain (e.g., < it >= tyr—), and u* = it*- < pu* > / < 
p >. (Throughout this work, volume averages are computed by multidimensional trapezoid 
rule, which is spectrally accurate whenever the fields are spatially periodic). The flow is 
then allowed to evolve without further forcing. The temporal evolution of total kinetic 
energy TKE = 0.5 < ptxj e 'ujj' > and turbulent Mach number M t are shown in Figs. 4 and 5, 
respectively. Energy spectra E(k) at selected times are compared with the initial spectrum 


= 166. 

= 1.0 (23) 

= 300 K 
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in Fig. 3. Here, k = |k| = \Jk\ + k\ + k%, where k = [fci, & 2 , A^] is the wavenumber vector in 
Fourier space. 

6 DNS: A Priori Testing 


The governing equations (Eq. 1) are solved by efficient pseudospectral methods using the 
algorithm developed by Obrist [13] for his Diploma thesis at ETH under the direction of 
the second author. The uniformly-spaced computational grid uses 96 points (48 Fourier 
modes) in each coordinate dimension. Because the flow is well resolved at the 96 3 resolution 
for the particular parameter values of interest, no explicit dealiasing is necessary (nor is it 
desirable). Time is advanced fully explicitly by the third-order low-storage Runge-Kutta 
method of Williamson [14]. 


As noted in Blaisdell et al. [12], the Reynolds number of the simulation is insufficiently 
high to result in a distinct inertial subrange with a characteristic -5/3 power-law decay. The 
temporal evolution of the spectrum of the present simulation is shown in Fig. 3. All analyses 
are performed at dimensionless time t = 1.51. For completeness, Figs. 6, 7, and 8 show the 
temporal evolution of the turbulent Reynolds number Rer, the Taylor microscale An, and 
the Taylor microscale Reynolds number Re\, which are defined, respectively, as follows: 


Rex — R& 


P = 


€ 


< > 2 

<p}j> 

<p> 

1 du'L 
Te < ° ii dx l > 


(25) 


A 2 


kl 


< > 

< (& 2 > 


(26) 


Re\ — Re < p> qXnji 
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(27) 



(M) 

A« 

(1,1) 

0.3579 

(1,2) 

0.2925 

(1,3) 

0.2949 

w 

0.2832 

w 

0.3600 

(2.3) 

0.2928 

n 


0.2882 

1 3?T 

0.2922 

(3,3) 

0.3566 


Table 1: Taylor-microscale components at dimensionless time t = 1.51. 



< Pu'M > 


T l/2 


< p > 


Note that q is a turbulent velocity scale, for which TKE = q 2 / 2, and that c is the turbulent 
dissipation rate. Like Blaisdell et al. [12] we observe two distinct microscales. The diagonal 
components of A# are consistent as are the off-diagonal components, as shown in Table 1. In 
general, the diagonal components are larger by 20-40 percent. Also as observed by Blaisdell 
et al. [12], the microscale diminishes at first and then grows gradually. For comparisons with 
the results of Blaisdell et al. [12], for which time is typically given in units U$/{TKE)q> we 
note that our dimensionless time should be scaled by a factor of approximately two. 


Fig. 9 presents contours of constant r u in an arbitrarily selected x-y plane that corre- 
sponds to z = 0.6545. The residual-stress component is computed exactly by grid-filtering 
the DNS results according to Eq. 11. Fig. 9 further compares the exact residual stress com- 
ponent with contours of the corresponding components of the Smagorinsky (Ml, Eq. 15), 
similarity (M2, Eq. 18), and generalized Clark models (M3, Eq. 22), respectively. For M2, 
which requires twice filtered quantities, the cutoff values for the grid and test filters are 
ot( — 7 t /2 and a c — 0.93, respectively. Fig. 9 clearly reveals model M2 to be much more 
highly correlated to the exact residual stress (E) than either Ml or M3. 


The numerical values of these correlations are summarized in Table 2 for all components 
of Tjy. In particular, Table 2 presents correlations of the form C(E, M) between the exact 
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(1-91, tt/2) 

C(E,M 1) 

C(E,M 2) 


ni 

0.066 

0.977 

0.619 

Tl2 

0.182 

0.975 

0.650 

Tl3 

0.125 

0.978 

0.546 

T22 

0.081 

0.976 

0.609 

T23 

0.113 

0.978 

0.540 

T 33 

0.214 

0.980 

0.412 

(tt/2,0.93) 

C(E,M 1) 

C(E,M 2) 

C(E,M 3) 

Til 

0.067 

0.910 

0.474 

Tl2 

0.183 

0.906 

0.518 

Tl3 

0.116 

0.913 

0.432 

T22 

0.082 

0.908 

0.461 

T23 

0.102 

0.913 

0.427 

T 33 

0.199 

0.921 

0.314 


Table 2: Correlations between exact (E) residual-stress fields and their modeled ( M ) 
counterparts, for case of differing grid and test filters, whose respective parameters ot c are 
provided in upper left box of each block. 


( E ) and modeled ( M ) residual-stress fields for the three models under consideration. Here, 
we define the correlation coefficient C between any two fields u and v in the customary way 
as follows: 


n( s < uv > <uXu> . . 

C(u ’ “ [(< u 2 > - < u > 2 )(< v 2 > - < t; > 2 )] 1 / 2 ( } 

In general, we observe poor correlations (typically C < 0.2) for Ml, moderately high corre- 
lations ( C « 0.6) for M3, and remarkably high correlations (C > 0.9) for M3. Two filtering 
operations are required to evaluate M2. Thus, the two values that appear in the upper left 
box of each block of Table 2 correspond to the respective cutoff values for the grid and test 
filters. 


Fig. 10 and Table 3 present similar results for the residual-stress components Qk (Eq. 14) 
with their modeled (Ml, M2) counterparts MQk (Eq. 17) and q * (Eq. 21), respectively. The 
same levels of correlations observed for components of the SGS-stress tensor are observed 
also for the thermal stresses. 


It is well known that scalar eddy viscosity models (e.g., Ml) do not perform well in 
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(1.91,i/2) 

C(E,M1) 

C(E, M2) 

Q i 

0.106 

0.977 

Q 2 

0.109 


Q 3 

0.236 

0.974 

(tt/2,0.93) 

C(E, Ml) 

C(E,M2) 

Q 1 

0.104 

■KEEfl 

Q 2 

0.110 

0.911 

Qz 

0.213 

0.912 


Table 3: Correlations between exact (E) energy residual fields Qk and their modeled ( M ) 
counterparts, for case of differing grid and test filters, whose respective parameters a c are 
provided in upper left box of each block. 

anisotropic flows (e.g., Compton and Eaton [15]) because the shear stress and the mean 
strain are misaligned. However, the component- wise variation of the correlations C(E, Ml) 
in Tables 2 and 3 comes initially as a surprise. For isotropic turbulence, one expects the 
correlations to be nearly the same across all components. Statistical isotropy, however, does 
not im ply instantaneous isotropy. Numerical experimentation with different random initial 
conditions and a priori analyses at different snapshots in time reveal these variations to be 
manifestations of instantaneous anisotropy. On the other hand, the similarity model (M2) 
appears relatively insensitive to anisotropy. 

Liu et al. [4] found Bardina’s [10] model to yield the highest correlation coefficients. 
Accordingly, Fig. 11 presents contours of the exact (E) and similarity (M2) modeled (2,3) 
components of the residual-stress tensor for the case in which the grid and test filters are 
identical, each with a c = 7t/2. In this case, the correlation C(E, M2) attains a surprisingly 
high value, in excess of 0.99. Table 4 below summarizes the correlation data for various values 
of the cutoff a c , which is given in the upper left box of each block in the table. The correlation 
coefficients of the form C(E, M2) range from a low of 0.86, for an extremely dissipative filter, 
to 0.999, for a very mild filter. As expected, the less dissipative the filter (i.e., the higher 
k c ), the higher the correlation coefficient. In addition to correlation coefficients, Table 4 
also indicates the ranges (maxima and minima) of the exact and modeled residual stresses. 
From these extrema, we make two important observations. First, the magnitudes of the 
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components of the exact residual-stress tensor r w depend significantly on the value of the 
cutoff a c . Second, whenever identical grid and test filters are exploited, r kl is apparently well 
approximated by Eq. 18 with cl « 1.0, provided the filter is hot too dissipative. Tables 2 
and 4 reveal that the fidelity of all three models appears to deteriorate as the filter becomes 
highly dissipative (although the correlations for Ml are so poor that any effect due to filter 
tuning is difficult to assess). 

To provide additional insight, Table 5 presents correlations of the form C(E, M) for the 
case where the grid filter is of fourth order, as discussed in Section 2. For M2 the test filter is 
of second order, so that the grid and test filters are of differing orders. The values shown in the 
upper left box of each block of Table 5 correspond to the cutoff values a e for the grid and test 
filters, respectively. In general, component- wise correlations for the generalized Smagorinsky 
model appear to drop below their previously low values, as suggested by a comparison of 
Tables 2 and 5. Moreover, correlations of the type C(E, M2) drop dramatically to range 
from 0.5 to 0.7 depending upon the values of the grid- and, test-filter cutoffs. The poorest 
correlations for M2 arise when the respective cutoffs are different. Curiously, M3 correlates 
as well or better for a fourth-order grid filter as it did for a second-order grid filter. This 
may suggest that the fourth-order term in Eq. 22 is significant. 

Finally, we consider the effect of the order property of the grid filter on the exact residual- 
stress tensor. Figure 12 compares the exact (2,2) component of t*i for grid filters of different 
orders and cutoffs. The first grid filter is of second order with a c = 1.16; the second 
is of fourth order with a c — 1.71. Contours are displayed at z = 0.6545. Clearly, the 
stress distribution is grid- filter dependent. To quantify this dependence, Table 6 presents 
correlation coefficients of the form C(E2,i?4), where £2 and E\ refer to exact quantities 
computed by second- and fourth-order grid filtering, respectively. All components show the 
same level of correlation, nominally 0.54. We conclude that the exact residual stresses depend 
on the selection of the grid filter in both their distribution and their scale. 
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C(£,M2) 

-Emax 


M 2 ma x 

■H 

I'll 

0.997 

0.112E-02 

-.654E-02 

0.896E-03 

-.572E-02 

r 12 

0.997 

0.315E-02 

-.288E-02 

0.271 E-02 

-.247E-02 

Tl3 

0.997 

liiraSFi 

-.291E-02 

0.298E-02 

-.246E-02 

T 22 

0.997 

0.954E-03 

-.728E-02 

0.782E-03 

-.614E-02 

^23 

0.997 

0.270E-02 


0.240E-02 

-.200E-02 

7 33 

0.997 

0.613E-03 


0.433E-03 

-.489E-02 

a c = 7r/2 

C(£,M2) 

■■FOMi 


M2max 


Til 

0.992 

0.00 


0.00 

-.125E-01 

Tl2 

0.992 

0.772E-02 

w 

0.620E-02 

-.512E-02 

Tl3 

0.993 

0.847E-02 

-821E-02 

0.653E-02 

-.469E-02 

T 22 

0.992 

0.00 

-.167E-01 

0.00 

-.122E-01 

T23 

0.993 

0.664E-02 

BSKW>B0ti 

0.516E-02 

-.432E-02 

T 33 

0.994 

0.00 

>.138E-01 

0.00 

! -.107E-01 

a c = 0.93 




M2max 

■Bl 

Til 

0.970 

HKEi 

-.315E-01 

0.00 

-.188E-01 

Tl2 

0.969 



0.890E-02 

-.854E-02 

Tl3 

0.975 

0.168E-01 


0.105E-01 

-.666E-02 

T22 

0.968 

0.00 


0.00 

-.210E-01 

T23 

0.975 

0.133E-01 

N4JOIE-OI 

0.900E-02 

-.624E-02 

T 33 

0.978 

0.00 

-.270E-01 

0.00 

-.162E-01 

a c — 0.45 


^max 

m 

M2max 


Til 

0.861 


rtEEKiaiCT 

HHESS1 

-.186E-01 

T 12 

0.872 

mi 

ISM12M! 

msrnmM 

-.911E-02 

Tl3 

0.887 

■iMiwaam 


um&m 


T 22 

0.865 


mm a 

0.00 

-.285E-01 

T23 

0.890 

0.217E-01f 

W..143E-01 

0.974E-02 

-.888E-02 

T 33 

0.901 



0.00 

-.226E-01 


Table 4: Correlations between exact (E) z$$4tyal fields and their similarity modeled (M2) 
counterparts, and maximum and minimum values of E and M2 fields, for case of identical 
grid and test filters, whose parameter a c is provided in upper left box of each block. 
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■bki 
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T 33 


ranr.rzi 

mbeeebi 

gflFMiiai 

11 

EH^®I 

j— 

7i2 

WBES^ 



ns 

buksii^ 


0.550 

. : ,i 



mebhi 

^23 

— m 

0.587 

—™™ 

ns lL. 

■linn 

i^essi 

liwfm 



WMM 

isarai 

ni 


iffipf 


n 2 

IK2U1 

■KM^I 

HHE222II 

na 

■HB21 


warn wt 

n 2 

■awrnnygi 

■K^l 

0.626 

7 23 


■K2^] 

BHBUiJMifl 

ns 

0.176 

mm 

0.425 


Table 5: Correlations between exact (E 1 ) residualfields and their modeled (Af) counterparts, 
for case with fourth-order grid filter and second-order test filter, whose respective cutoffs a c 
are provided in upper left box of each block. 
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(1.16,1.71) 

C(E2, E4) 

Til 

0.548 

T"l2 

0.548 

Tl3 

0.536 

7*22 

0.546 

r 23 

0.538 

T33 

0.544 

Qi 

0.537 

Q 2 

0.539 

Q 3 

0.535 


Table 6: Correlations between exact residual-stress components computed by second-order 
(E2) and fourth-order (E4) grid filtering, whose respective cutoffs a c are given at upper left. 

7 Discussion 


In the conventional practice of LES, the grid filter and the SGS-model have often been 
treated as if they were independent choices (as observed by Piomelli et al. [16]). In contrast, 
a priori analyses with a continuously tunable filter reveal that the filter, the exact SGS 
stresses, and the model are closely interrelated. Thus, the filter and the model should not 
be specified independently. 

For reasons of practicality, our computational study was limited to a single test flow at 
a relatively low Re\, much lower, for example, than the value Re\ = 310 in the experiment 
of Liu et al. [4]. It is natural to wonder if the high correlations observed for the stress- 
similarity model (M2) generalize to other flows and to higher Re\. Several considerations 
suggest an affirmative answer. First, correlations of 0.85 were observed for M2 in the phys- 
ical experiment of Liu et al. [4] whenever identical second-order grid and test filters were 
exploited. Second, by comparing Fig. 8 and Table 7 below, one observes some dependence of 
the correlation coefficient on Re a, but it appears to be relatively weak. Third, with identical 
second-order Pade grid and test filters each of ot c = tt/ 2, Adams et al. [17] have recently 
observed correlations well in excess of 0.90 for DNS of a forced turbulent Mach 3.0 boundary- 
layer flow along a compression ramp, for Which Re\ grows from a value of 47 in the incoming 
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t 

C{E,M 2) 

.127 

0.966 

.287 

0.971 

.644 

0.978 

.829 

0.982 

.977 

0.985 

1.08 

0.987 

1.18 

0.988 

1.24 

0.989 

1.40 

0.991 

1.51 

0.992 


Table 7: Time evolution of correlation coefficients of the form C(r n> £ u ) for identical grid 
and test filters with a c = 7r/2. 


boundary layer to 326 midway up the compression ramp. Moreover, even with a much more 
dissipative filter (a c = 0.53), correlations on the order 0.8 are observed. It is important to 
note that, in contrast to the present test case, the flow considered by Adams et al. is highly 
anisotropic, highly compressible, and wall-bounded. 


8 Conclusions 

The decay of isotropic turbulence in a compressible flow is examined by direct numerical 
simulation (DNS). A priori analyses of the DNS data are then performed to evaluate three 
candidate subgrid-scale stress (SGS) models for large-eddy simulations (LES): the generalized 
(compressible) Smagorinsky model (Ml) of Erlebacher et al. [8], the stress-similarity model 
(M2) of Liu et al. (4], and the generalized Clark model (M3) of Vreman et al. [7]. 

The focus of the study is an examination of the effect of the grid and test filters on the 
fidelity of the models. For this purpose, tunable, one-parameter filters of Pade type are 
exploited. Tunable filters are desirable for LES in that they allow the cutoff wavenumber 
k c to be specified independently of the grid resolution Ax. In our view, the cutoff should 
be selected on the basis of physical considerations and the grid resolution on the basis of 
mathematical ones. 
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In the a priori analyses, the modeled residual stresses are Correlated against their exact 
counterparts for selected values of the cutoff wavenumber k c of the filter. On the basis of 
these results, the following conclusions are offered: 

1. The distribution and the scale of the exact residual-stress tensor r w depends strongly on 

both the cutoff k c and the order property of the grid filter. This dependence is partially 
revealed by the simple recognition that the residual stresses must vanish entirely in the 
limit as k c -»• k v , where k v is the Kolomogorov wavenumber. Consequently, in general, 
no SGS-stress model can be independent of the grid filter. In particular, there is no 
universal Smagorinsky constant; the constant necessarily depends of the character of 
the grid filter. n;^ 

2. Correlations of the form C(E , M) between the exact (E) and modeled (M) SGS stresses 
are typically very high ( C > 0.9) for the stress-similarity model (M2), moderately good 
(C « 0.6) for the Clark model (M3), and poor ( C < 0.2) for the Smagorinsky model 
(Ml). 

3. For each of the three models, higher correlations C(E,M) are observed for less dissi- 
pative grid filters and vice versa. This is an expected result because (for a fixed grid 
resolution) dissipative filters are characterized by relatively low cutoff wavenumbers k c , 
which in turn results in more energetic subgrid scales of motion. 

4. Model M2 performs remarkably well in a priori tests provided the grid and test filters 
have the same order property. The highest correlations, well in excess of 0.90, are 
obtained for M2 whenever the grid and test filters are identical. If the grid and test 
filters are identical and the filters are not too dissipative, model M2 appears to capture 
not only the correct distribution of stress but also the correct amplitude whenever the 
model constant Ci = 1.0. This result numerically replicates Bardina’s [10] result, albeit 
for compressible flow. 

5. The present results for models Ml and M2 are consistent with the experimental find- 
ings of Liu et al. [4], who did not consider M3. 
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6. Whereas models Ml and M3 appear to degrade whenever the flow is instantaneously 
anisotropic, model M2 appears relatively insensitive to anisotropy. 

As a concluding thought, we acknowledge the well-known experience that both Bardina’s 
model (e.g., Liu et al. [4]) and Clark’s model (Vreman et al. [5]) are insufficiently dissipative 
for practical applications to LES. However, in our view, the physics and the mathematics 
of LES need decoupling. It appears increasing likely that stress-similarity models (M2, for 
example) capture the physics (e.g., backscatter) with high fidelity. That such models “blow 
up” in practice is not surprising in light of Fig. 1. The physical viscosity, which dampens 
the highest wavenumbers in DNS, is ineffective in LES at the cutoff wavenumber k c . We 
conclude that the instability of similarity models is ultimately a mathematical issue, not a 
physical one. 
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Appendix: One-Parameter Pade Filters 


Let f be a vector containing n + 1 values fi = f{xi) obtained by sampling a continuous 
function f(x) at a set of equally spaced nodes x, = xo + iAx, i = 0, 1, n. Further let 
f = Lf denote the vector of filtered values obtained by applying the linear filter L to f. 
Here, we consider Pade filters, for which L = M~ l E , whereby Mi = Ef, and M and E are 
banded matrices. 


Second-Order Pade Filters 


Following Lele [3], a second-order, one-parameter (C) Pade filter is constructed by considering 
the symmetric pointwise scheme 

C/*-l + /* + C/i+l = a fi + 2 ^* -1 "** ^ <+1 ^ ’ (* = — n — 1) (29) 

Various treatments are possible for the boundary nodes i = 0 and i = n. However, the sim- 
plest approach, sufficient for our present purposes, is to impose no filtering at the boundaries. 
For the second-order Pade scheme, M and E are tridiagonal matrices. 


It is instructive consider the action of the filter of Eq. 29 on the single Fourier mode 
exp (tkx) (i = y/—l), from which we obtain the complex transfer function 

H(a) = a + bcos (?) (30) 

w 1 + 2£ cos(a) 

where a = fcArr. For a fixed grid increment Ax, the transfer function defines the action of 
the filter on modes of wavenumber k. For applications to LES, we consider only low-pass 


filters, for which H(0) = 1 and H( n) = 0. In combination, these two constraints imply that 
b = a = 0.5 + C- Admissible values of the parameter are — | < C < Whenever C = 0> 
M is the identity matrix, and the filter is said to be fully explicit. The fully explicit case 
corresponds to a discrete top-hat filter with weights at adjacent nodes of [1/4, 1/2, 1/4]. 
Otherwise, the filter is implicit. The value C = 2 yi e ^ s M = E, which turns off the filter. 
For all admissible parameter values, the matrices M and E are diagonally dominant with 
positive diagonal elements, in which case the operator L is positive semidefinite. 
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By analogy to discrete differentiation' operators, to each discrete filter is associated a 
Taylor-series expansion of the form 

f(xi) = f(ii) + ai Ax f'(xi) + a 2 Ax 2 f"(xi) + a 3 A x z f"'(xi) + ... (31) 

Here, for notational simplicity, the factorials in the Taylor expansion have been absorbed 
into the coefficients Oj. By applying the Taylor-series representation of the filter to exp(ikxi), 
we obtain the corresponding Taylor series of the transfer function, namely 

H(a) = 1 + aim - o 2 a 2 - a 3 wJ» 3 •+■ ... (32) 


Equating successive derivatives of the two forms of the transfer function (Eqs. 30 and 32) at 
a = 0, yields the coefficients Cj of the transfer function. Thus, the Taylor-series expansion 
of the filter implies the transfer function and vice versa. For the Pade filter of interest, the 
first two coefficients are 


a x = H\ 0) = 0 

-2* = >(0)='-i + r ^ (33) 


In general, a filter is said to be of order m if the first non-vanishing coefficient of its 
Taylor-series expansion is a*, . Thus, the Fade scheme is of second order (provided ( # 0.5). 
Equivalently, a filter is of order m if all derivatives of the transfer function of order less than 
m vanish at the origin. In general, filters associated with symmetric stencils (and coefficients) 
are of even order with purely real transfer functions. As a consequence, they filter without 
phase error. 


As clearly indicated by Eq. 33, the coefficients of the Taylor-series expansions Eq. 31 and 
Eq. 32 are functions of the parameter. To quantify this dependence, we define at c , the cutoff 
wavenumber, such that 


, C f= H(a c ) = 


(l + 2C)(l+cosa e ) 


(34) 


2(l + 2<cosa c ) 

where 0 < C < 1 but is otherwise arbitrary. (In signal processing, it is customary to define 
the cutoff such that 1/2 = |i?(a^)| 2 .) Eq. 34 relates the original parameter C to the cutoff 
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wavenumber 0 < a c < n. The (unconventional) choice C = i/2 yields the particularly useful 
relationship 2£ = -cos(a c ), which in turn yields 

(1 — cos a c )(l + cos a) 


H{a , a c ) = 


(35) 


2(1 - cos a c cos a) 

Fig. 2 in Section 3 compares the transfer functions of the variable Pade filter for selected 
values of the cutoff a c . 


Fourth- Order Pade Filters 


We adapt the three-parameter, symmetric, fourth-order filter of Lele [3] whose M and E 
matrices are penta- and hepta-diagonal, respectively. The pointwise scheme for nodes * = 
2,3, ...,n - 2 is 

0fi-2+Cfi-\ + fi+Cfi+l+Pfi+2+ = a fi + 2 (^‘ _1 +/*+l) + 2 + 2 ^‘“ 3 ^* +3 ^ 

where 

a = +(5 + 6C-6/?+16d)/8 

b * +(1 + 2C + 20 - 2d)/2 (37) 

c = -(1 - 2C - 14/3 + 16d)/8 

(38) 


By enforcing c — d = 0, one arrives at a one-parameter, fourth-order filter whose matrices 
M and E are penta-diagonal and tridiagonal, respectively, and for which 0 = (1 - 2£)/14. 
The stencil is closed by imposing no filtering at the boundaries and by using the explicit 
fourth-order formula below at the gridpoint adjacent to the boundary point. 

/i = ^/i + jg(/o + 6/2 — 4/3-I-/4) (39) 

An analogous formula is used at gridpoint n — 1. Admissible values of the parameter are 
— | < C < (Gerschgorin’s theorem guarantees positivity only for — | < C < §•) A 
linear relationship between the original parameter C and — 1 < cos(a c ) < 1 establishes 
12C = -[1 + 7cos(a c )], which defines the cutoff wavenumber a c in terms of C- 
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Figure 1: Kolmogorov energy spectrum for idealized LGS. 



Figure 2: Transfer function of one-parameter family of second-order low-pass filters of Pade 
type for selected values of dimensionless cutoff wavenumber a c . 
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Figure 3: Kinetic energy spectrum at selected times 



time, t 

Figure 4: Temporal evolution of total kinetic energy TKE. 
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Figure 7: Temporal evolution of Taylor microscale An. 
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Figure 8: Temporal evolution of Tayldr-microscale Reynolds number, Re\ 
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Figure 9: Contours of exact (E) and modeled (M) r n component of residual stress in 
selected x — y plane at t = 1.5, for grid- and test-filter cutoff values a c = 7r/2 and a c = 0.93, 
respectively. 
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Figure 10: Contours of exact (E) and modeled (M) qi component of residual thermal stress 
in selected x - y plane at t = 1.5, for grid-filter and test-filter cutoff values a c = 7r/2 and 
a c = 0.93, respectively. 
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Figure 11; Contours of exact (E) and modeled(M) r 23 component of residual stress in 
selected x - y plane at t = 1.5, for identical grid and test filters of cutoff values a c = ir/2. 
For ranges of E and M2 values, see Table 4. 
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Figure 12: Comparison of exact ( E ) residual-stress component r » for grid filters of second 
(E2) and fourth (E4) order. 
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